%% fig5a_modified.m
% 改进版 Fig.5：同时绘制分岔图和最大李雅普诺夫指数图，
% 其中李雅普诺夫指数的计算调用封装的 LE.m 函数。

clear; clc; close all;

%% 参数设置
b     = 0.5;
alpha = 2.1;      
beta  = 2;
gamma = 2;
Im    = 0.15;        % 外激励幅值
F     = 0.15;        % 外激励频率

% 控制参数 a 的取值范围
a_vals = linspace(0.57, 0.87, 200);

% 时间参数设置
T_final = 1200;    % 总积分时间
T_trans = 200;     % 舍去的瞬态时间
% 整个积分区间用于分岔图部分：用 [T_trans, T_final] 数据提取局部峰值
optsODE = odeset('RelTol',1e-9,'AbsTol',1e-9);

% 分段积分参数用于李雅普诺夫指数计算
segment = 5;      % 每段积分时间
% 注意：LE.m 内部会根据 T_trans, T_final, segment 计算 N_seg

% 预分配数据
all_a_bif    = [];  % 存储每个 a 对应的分岔图数据（重复 a）
all_Vx_peaks = [];  % 存储对应的局部峰值 Vₓ
lyap_vals    = zeros(length(a_vals),1);  % 每个 a 对应的最大李雅普诺夫指数

%% 针对每个 a 值进行计算（并行循环）
parfor i = 1:length(a_vals)
    a = a_vals(i);
    
    %% 1. 去除瞬态积分（计算扩展系统初值，用于后续两部分）
    % 初始状态：主系统 Y0=[q; sigma; phi]=[0;0;0]，扰动向量选取非零向量归一化
    Y0 = [0; 0; 0];
    v0 = [1; 0; 0]; v0 = v0/norm(v0);
    X0 = [Y0; v0];
    
    tspan_trans = [0, T_trans];
    [~, X_trans] = ode45(@(t,X) cnn(t, X, a, b, alpha, beta, gamma, Im, F), tspan_trans, X0, optsODE);
    X_current = X_trans(end,:)';  % 用于后续积分
    
    %% 2. 分岔图数据提取：积分区间 [T_trans, T_final]
    tspan_bif = [T_trans, T_final];
    [t_bif, X_bif] = ode45(@(t,X) cnn(t, X, a, b, alpha, beta, gamma, Im, F), tspan_bif, X_current, optsODE);
    % 提取主系统状态 Y=[q; sigma; phi]
    Y_bif = X_bif(:,1:3);
    q     = Y_bif(:,1);
    sigma = Y_bif(:,2);
    % 计算 Vₓ = (alpha+beta*sigma).*q
    Vx = (alpha + beta.*sigma).*q;
    % 利用 findpeaks 从稳态段中提取局部峰值（Vₓ 的局部最大值）
    [pks, ~] = findpeaks(Vx, t_bif);
    
    % 将该 a 下的分岔图数据保存
    a_temp = a * ones(length(pks),1);
    
    %% 3. 调用封装的 LE.m 函数计算最大李雅普诺夫指数
    % LE.m 的调用语法：
    % LE_val = LE(a, b, alpha, beta, gamma, Im, F, T_trans, T_final, segment)
    lyap_val = LE(a, b, alpha, beta, gamma, Im, F, T_trans, T_final, segment);
    
    %% 将数据赋值到对应输出变量
    lyap_vals(i) = lyap_val;
    all_a_bif = [all_a_bif; a_temp];
    all_Vx_peaks = [all_Vx_peaks; pks];
end

%% 绘图：在同一 figure 内绘制分岔图与李雅普诺夫指数图
figure;

% 子图1：分岔图（散点图）
subplot(2,1,1);
scatter(all_a_bif, all_Vx_peaks, 10, 'filled');
xlabel('控制参数 a');
ylabel('V_x 的局部最大值');
title('分岔图');
grid on;

% 子图2：最大李雅普诺夫指数曲线
subplot(2,1,2);
plot(a_vals, lyap_vals, 'b.-', 'LineWidth', 1.5);
xlabel('控制参数 a');
ylabel('最大李雅普诺夫指数');
title('李雅普诺夫指数');
grid on;
